#### "Gender Inequality and Authoritarian Regimes" ####
# authors: "Lars Pelke"
# date: 2020-01-20

## Set WD ##

## use 01_data_wrangling before running this file  OR

## load data ##

vdem <- readRDS("calculations/data/vdem_merged.rds")



#### Figures A1 Main Paper ####

vdem$gwf_type <- as.factor(vdem$gwf_type)
vdem <- vdem %>%
  drop_na(gwf_type)

ggplot(data = vdem , aes(x = gwf_type, y =v2xpe_exlgender)) +
  geom_violin(alpha=1, position = position_dodge(width = .75), size=1, color="black", fill = "grey90") +
  geom_boxplot(notch = TRUE, outlier.size = -1, color="black",lwd=1.2, alpha = 0.4, width = .2)+
  geom_jitter(shape = 16, size=1,  color="grey", alpha= 1) +
  theme_pubr() +
  labs(x = "Regime Types", 
       y = "Political Exclusion by Gender")

ggsave("calculations/Output/Figure1.pdf", height = 14, width = 24, units= c("cm"))


#### Figure 2 Main Paper ####

ggplot(data = vdem , aes(x = gwf_type, y =v2xps_party)) +
  geom_violin(alpha=1, position = position_dodge(width = .75), size=1, color="black", fill = "grey90") +
  geom_boxplot(notch = TRUE, outlier.size = -1, color="black",lwd=1.2, alpha = 0.4, width = .2)+
  geom_jitter(shape = 16, size=1,  color="grey", alpha= 1) +
  theme_pubr() +
  labs(x = "Regime Types", 
       y = "Party Institutionalization")

ggsave("calculations/Output/Figure2.pdf", height = 14, width = 24, units= c("cm"))

#### Prepare Data for analysis ####

summary(vdem$gender_inclusion_input_mean)
vdem$gender_inclusion_input_mean <- rescale(vdem$gender_inclusion_input_mean, to = c(0,1))
summary(vdem$gender_inclusion_input_mean)

summary(vdem$gender_inclusion_output_mean)
vdem$gender_inclusion_output_mean <- rescale(vdem$gender_inclusion_output_mean, to = c(0,1))
summary(vdem$gender_inclusion_output_mean)

vdem <- vdem %>%
  drop_na(cown, year, gwf_type, gwf_type_mp)

vdem2 <- vdem %>%
  drop_na(cown, year, gwf_type, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio)

vdem3 <- vdem %>%
  drop_na(cown, year, gwf_type, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio, v2xps_party)


#### Descriptive Tables ####

#functions 
mean.new <- function(v) {
  if (all(is.na(v))) { return(NA) } else { return(mean(v, na.rm=T)) }
}

## Table 2 Main Paper ##
table1 <- vdem %>%
  group_by(gwf_type_mp) %>%
  summarize(mean_input = mean.new(gender_inclusion_input_mean), 
            mean_output = mean.new(gender_inclusion_output_mean))%>%
  rename("Average WPI on\nthe input side" = mean_input, 
         "Average WPI on\nthe output side" = mean_output, 
         "Regime Type" = gwf_type_mp)

library(stargazer)
stargazer(table1, summary=FALSE)

## country examples for Table 1##

country_examples <- vdem %>%
  filter(gwf_type_mp=="personal regime with ME") %>%
  dplyr::select(country_name, year, gwf_type_mp, gwf_casename, gender_inclusion_input_mean, gender_inclusion_output_mean)

rm(table1, country_examples)

## Table 3 Main Paper ##

vdem <- vdem %>%
  ungroup() %>%
  mutate(gwf_type_mp_num = case_when(gwf_type_mp == "party regime with ME" ~ 1, 
                                     gwf_type_mp == "party regime" ~ 2,
                                     gwf_type_mp == "monarchy with ME"~ 3, 
                                     gwf_type_mp == "monarchy"~ 4, 
                                     gwf_type_mp == "personal regime with ME"~ 5, 
                                     gwf_type_mp == "personal regime"~ 6, 
                                     gwf_type_mp == "military regime with ME"~ 7, 
                                     gwf_type_mp == "military regime"~ 8,
                                     gwf_type_mp == "communist party with ME"~ 9,
                                     gwf_type_mp == "communist party"~ 10)) 
vdem <- vdem %>%
  group_by(country_name) %>%
  mutate(group_id = cumsum(gwf_type_mp_num != lag(gwf_type_mp_num, default = FALSE))) # group_id for each regime spell

vdem <- vdem %>%
  mutate(regime_type_id = str_c(country_name, group_id, sep = "_"))

table3 <- vdem %>%
  group_by(regime_type_id) %>%
  summarize(regime_type = first(gwf_type_mp), 
            start_year = first(year), 
            end_year = last(year), 
            mean_input = mean.new(gender_inclusion_input_mean))

table3a <- table3 %>%
  slice_max(order_by = mean_input, n = 10) %>%
  rename("Regime Type" = regime_type, 
         "Begin" = start_year, 
         "End" = end_year, 
         "Mean WPI on the input side" = mean_input)

stargazer(table3a, summary=FALSE)

table3b <- table3 %>%
  slice_min(order_by = mean_input, n = 10) %>%
  rename("Regime Type" = regime_type, 
         "Begin" = start_year, 
         "End" = end_year, 
         "Mean WPI on the input side" = mean_input)

stargazer(table3b, summary=FALSE)


## Table 4 Main Paper ##

table4 <- vdem %>%
  group_by(regime_type_id) %>%
  summarize(regime_type = first(gwf_type_mp), 
            start_year = first(year), 
            end_year = last(year), 
            mean_output = mean.new(gender_inclusion_output_mean))

table4a <- table4 %>%
  slice_max(order_by = mean_output, n = 10) %>%
  rename("Regime Type" = regime_type, 
         "Begin" = start_year, 
         "End" = end_year, 
         "Mean WPI on the output side" = mean_output)

stargazer(table4a, summary=FALSE)

table4b <- table4 %>%
  slice_min(order_by = mean_output, n = 10) %>%
  rename("Regime Type" = regime_type, 
         "Begin" = start_year, 
         "End" = end_year, 
         "Mean WPI on the output side" = mean_output)

stargazer(table4b, summary=FALSE)



#### Computing group_mean and de-meaned variables ####

vdem2 <- sjmisc::de_mean(vdem2, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, 
                         grp = cown)

vdem3 <- sjmisc::de_mean(vdem3, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, v2xps_party,  
                         grp = cown)


#### Model Building: Complex REWB Models ####

## Relevel Factor's ##
vdem$gwf_type <- as.factor(vdem$gwf_type)
vdem<- vdem %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))

vdem$gwf_type_mp <- as.factor(vdem$gwf_type_mp)
vdem <- vdem %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))

vdem2$gwf_type <- as.factor(vdem2$gwf_type)
vdem2 <- vdem2 %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))

vdem2$gwf_type_mp <- as.factor(vdem2$gwf_type_mp)
vdem2 <- vdem2 %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))

#### Model Building step by step ####
## compare https://strengejacke.github.io/mixed-models-snippets/random-effects-within-between-effects-model.html 

## Growth Curve model

growth_model <- lmer(v2xpe_exlgender ~ year + gwf_type + v2xps_party + 
                       (1 + year | country_name),
                     data = vdem3, 
                     REML = TRUE) # REML estimation for parameters
summary(growth_model)
isSingular(growth_model, tol = 1e-05) # FALSE

## Complex REWB ##

complexREWB<- lmer(v2xpe_exlgender ~ year + gwf_type + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name),
          data = vdem3, 
          REML = TRUE) # REML estimation for parameters
isSingular(complexREWB) # FALSE

complexREWB.1 <- lmer(v2xpe_exlgender ~ year + gwf_type + v2xps_party_dm + v2xps_party_gm + 
              (1 + year | country_name) + (0 + v2xps_party_dm | country_name),
            data = vdem3,  # assume independence between random slopes and no coveriance
            REML = TRUE) # REML estimation for parameters
isSingular(complexREWB.1, tol = 1e-05) # FALSE

## Simple REWB model ##

simpleREWB <- lmer(v2xpe_exlgender ~ year + gwf_type + v2xps_party_dm + v2xps_party_gm + 
                     (1 + year | country_name), 
                   data = vdem3,
                   REML = TRUE) # REML estimation for parameters
isSingular(simpleREWB, tol = 1e-05) # FALSE

## Mundlak model ##

mundlak <- lmer(v2xpe_exlgender ~ year + gwf_type + v2xps_party + v2xps_party_gm + 
                  (1 + year | country_name), 
                data = vdem3,
                REML = TRUE)
isSingular(mundlak, tol = 1e-05) # FALSE

#### Comparison of Models and check suitable model ####

tab_model(growth_model, complexREWB, complexREWB.1, simpleREWB, mundlak,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Classical RE", "Complex REWB", "Complex REWB (2)", "Simple REWB", "Mundlak"))

## REWB or simple RE-model ? ##
anova(growth_model, mundlak) # significant improvement of the Mundlak-model over the simple 
# RE-model, indicating that it makes sense to model within- and 
# between-subjects effects

## Simple or Complex REWB? ##
anova(simpleREWB, complexREWB)
# Complex REWB Model: significant improvement of the Complex REWB over the the simple REWB model. 

## REWB model or FE model ? ##

fixed <- felm(v2xpe_exlgender ~ year + gwf_type + v2xps_party_dm  | country_name,
              data = vdem3)

tab_model(fixed, growth_model, complexREWB,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("FE-model with ID", "Classical RE", "Complex REWB"))

## fixed or Complex REWB? ##

anova(fixed, complexREWB) # no significant improvement of the complex REWB over the simple 
# REWB-model, indicating that it makes no sense to model random slopes and covariance

#### Residual diagnostics for complex REWB model ####

## Package ## 
library(DHARMa)
citation("DHARMa")

## Calculating scaled residuals ##

simulationOutput <- simulateResiduals(fittedModel = complexREWB, n = 500)
simulationOutput <- recalculateResiduals(simulationOutput, group = vdem3$country_name)

## Plotting the scaled redisuals ##

plot(simulationOutput, quantreg = FALSE) # plot leads to the conclusion that there is no big problem with outliers or uniformity

## Goodness-of-fit tests ##

# Uniformity: ok
testUniformity(simulationOutput = simulationOutput)

# Over/underdispersion: no problem 
testDispersion(simulationOutput)

# Heteroscedasticity: no problem
plot(simulationOutput, quantreg = FALSE) 

rm(mundlak, simpleREWB, complexREWB, complexREWB.1, fixed, growth_model, curWarnings, simulationOutput)



############################################################################################################
############################################################################################################

#### Bayesian Factor Analysis Political Exclusion by Gender Index #####

library(MCMCpack)

## Bayesian Factor Analysis ##

posterior_vdem_factor <- MCMCfactanal(~v2pepwrgen + v2clgencl + v2peapsgen + v2peasjgen + v2peasbgen, factors=1,
                                   verbose=0, store.scores=FALSE, a0=1, b0=0.15,
                                   data=vdem, burnin=1000, mcmc=20000, thin=1)

## First, the loadings (Lambda). ##

loadings.pol_inst <- summary(posterior_vdem_factor)$statistics[1:5]
loadings.pol_inst
ci.loadings.pol_inst <- summary(posterior_vdem_factor)$quantiles[1:5,c(1,5)]
ci.loadings.pol_inst

## Second, Uniqueness ##
uniquenesses.pol_inst <- data.frame(summary(posterior_vdem_factor)$statistics[6:10])
names(uniquenesses.pol_inst)[1] <- "uniquenesses"
uniquenesses.pol_inst


#### Regime Types Changes within countries ####

vdem_reg_change <- vdem %>%
  group_by(country_name) %>%
  mutate(reg_change = ifelse(gwf_type != lag(gwf_type), 1, 0)) %>%
  dplyr::select(country_name, year, gwf_type, reg_change)

table(vdem_reg_change$reg_change)


#### Regime Types Changes within countries ####

vdem_reg_change <- vdem %>%
  group_by(country_name) %>%
  mutate(reg_change = ifelse(gwf_type_mp_num != lag(gwf_type_mp_num), 1, 0)) %>%
  dplyr::select(country_name, year, gwf_type, reg_change)

table(vdem_reg_change$reg_change)

############################################################################################################
############################################################################################################

#### Main Models, M1 - M4 ####

m1<- lmer(v2xpe_exlgender ~ year + gwf_type + 
            (1 + year | country_name), 
          data = vdem, 
          REML = TRUE)
sjstats::icc(m1)

m2 <- lmer(v2xpe_exlgender ~ year + gwf_type_mp + 
             (1 + year | country_name), 
           data = vdem, 
           REML = TRUE)

m3<- lmer(v2xpe_exlgender ~ year + gwf_type + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1  + year | country_name), 
          data = vdem2, 
          REML = TRUE)

m4<- lmer(v2xpe_exlgender ~ year + gwf_type_mp + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1 + year | country_name), 
          data = vdem2, 
          REML = TRUE)

#### Extracting findings and plotting effects ####

tab_model(m1, m2, m3, m4,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3", "Model 4"))

library(texreg)
texreg(list(m1, m2, m3, m4), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Communist Party Regime",
                             "Monarchy",
                             "Party Regime",
                             "Personal",
                             "Communist Regime",
                             "Communist Regime with MC", 
                             "Military Regime with MC", 
                             "Monarchy", 
                             "Monarchy with MC", 
                             "Party Regime", 
                             "Party Regime with MC", 
                             "Personal Regime", 
                             "Personal Regime with MC", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Womens's Political Exlcusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "tab:tab-01", 
       longtable = TRUE)

#### Dot-Whisker Plot ####

library(dotwhisker)
library(broom.mixed)


m1_df <- tidy(m1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m2_df <- tidy(m2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m3_df <- tidy(m3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
             conf.method = "Wald")
m4_df <- tidy(m4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
             conf.method = "Wald")

m1_df <- m1_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal")) %>%
  filter(term!= "Intercept" & term!="Year")

m2_df <- m2_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC")) %>%
  filter(term!= "Intercept" & term!="Year")

m3_df <- m3_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal", 
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Monarchy" | term == "Party" | term =="Personal")

m4_df <- m4_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC", 
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Communist Party with MC" | term == "Military Regime with MC" | 
           term =="Monarchy" | term == "Monarchy with MC" | term=="Party Regime" | term == "Party Regime with MC" | 
           term =="Personal Regime" | term == "Personal Regime with MC")

m_models <- rbind(m1_df, m3_df)
plot_dot_whisker <- {dwplot(m_models, 
                          dot_args = list(size = 3, aes(shape = model)),
                          whisker_args = list(size = 1, aes(linetype = model)),
                          dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.1),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 1", "Model 3"),
                      labels = c("Model 1", "Model 3")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 1", "Model 3"),
                         labels = c("Model 1", "Model 3"))
  } 
ggsave("calculations/Output/Figure3.pdf", height = 14, width = 24, units= c("cm"))

m_models <- rbind(m2_df, m4_df)
plot_dot_whisker <- {dwplot(m_models,
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr()  + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.85, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 2", "Model 4"),
                      labels = c("Model 2", "Model 4")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 2", "Model 4"),
                         labels = c("Model 2", "Model 4"))
} 
ggsave("calculations/Output/Figure4.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m1)
car::vif(m2)
car::vif(m3)
car::vif(m4)



####  Models with Party Institutionalization as Indepedent Variable ####

## Relevel Factors ##
vdem3$gwf_type <- as.factor(vdem3$gwf_type)
vdem3 <- vdem3 %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))

vdem3$gwf_type_mp <- as.factor(vdem3$gwf_type_mp)
vdem3 <- vdem3 %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))


m5<- lmer(v2xpe_exlgender ~ year + gwf_type + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)
isSingular(m5, tol = 1e-05)

m6<- lmer(v2xpe_exlgender ~ year + gwf_type_mp + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)

m7 <- lmer(v2xpe_exlgender ~ year + gwf_type + v2xps_party_dm + v2xps_party_gm + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
             e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
             e_total_resources_income_pc_gm + 
             (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
           data = vdem3, 
           REML = TRUE)

m8<- lmer(v2xpe_exlgender ~ year + gwf_type_mp + v2xps_party_dm + v2xps_party_gm + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)

#### Texreg Findings ####

tab_model(m5, m6, m7, m8, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 5", "Model 6", "Model 7", "Model 8"))

texreg(list(m5, m6, m7, m8), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Communist Party Regime",
                             "Monarchy",
                             "Party Regime",
                             "Personal",
                             "Party Institutionalization", 
                             "Party Institutionalization (b)", 
                             "Communist Regime",
                             "Communist Regime with MC", 
                             "Military Regime with MC", 
                             "Monarchy", 
                             "Monarchy with MC", 
                             "Party Regime", 
                             "Party Regime with MC", 
                             "Personal Regime", 
                             "Personal Regime with MC", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Womens's Political Exlcusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "tab:tab-01", 
       longtable = TRUE)


#### Dot-Whisker Plots ####

m5_df <- tidy(m5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m6_df <- tidy(m6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m7_df <- tidy(m7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m8_df <- tidy(m8, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

m5_df <- m5_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)")) %>%
  filter(term!= "Intercept" & term!="Year")

m6_df <- m6_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)")) %>%
  filter(term!= "Intercept" & term!="Year")

m7_df <- m7_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Monarchy" | term == "Party" | term =="Personal" |
           term == "Party Institutionalization" | term == "Party Institutionalization (b)")

m8_df <- m8_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 8") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC",
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Communist Party with MC" | term == "Military Regime with MC" | 
           term =="Monarchy" | term == "Monarchy with MC" | term=="Party Regime" | term == "Party Regime with MC" | 
           term =="Personal Regime" | term == "Personal Regime with MC" | term == "Party Institutionalization" | 
           term == "Party Institutionalization (b)")

m_models <- rbind(m5_df, m7_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 5", "Model 7"),
                      labels = c("Model 5", "Model 7")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 5", "Model 7"),
                         labels = c("Model 5", "Model 7"))
} 
ggsave("calculations/Output/Figure5.pdf", height = 14, width = 24, units= c("cm"))

m_models <- rbind(m6_df, m8_df)
plot_dot_whisker <- {dwplot(m_models,
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr()  + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 6", "Model 8"),
                      labels = c("Model 6", "Model 8")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 6", "Model 8"),
                         labels = c("Model 6", "Model 8"))
} 
ggsave("calculations/Output/Figure6.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m5)
car::vif(m6)
car::vif(m7)
car::vif(m8)


############################################################################################################
############################################################################################################

#### Input Side ####

##load data ##

vdem <- readRDS("calculations/data/vdem_merged.rds")

#### Prepare Data for analysis ####
library(scales)

summary(vdem$gender_inclusion_input_mean)
vdem$gender_inclusion_input_mean <- rescale(vdem$gender_inclusion_input_mean, to = c(0,1))
summary(vdem$gender_inclusion_input_mean)

summary(vdem$gender_inclusion_output_mean)
vdem$gender_inclusion_output_mean <- rescale(vdem$gender_inclusion_output_mean, to = c(0,1))
summary(vdem$gender_inclusion_output_mean)

# delete those countries with less than 3 country-year observations to dela with singular fit 

vdem <- vdem %>%
  drop_na(cown, year, gwf_type, gender_inclusion_input_mean)

vdem <- vdem %>% # delete those countries with less than 3 country-year observations to dela with singular fit 
  dplyr::group_by(cown) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(cown)

vdem2 <- vdem %>%
  drop_na(cown, year, gwf_type, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio)

vdem3 <- vdem %>%
  drop_na(cown, year, gwf_type, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio, v2xps_party)

#### Computing group_mean and de-meaned variables ####

vdem2 <- sjmisc::de_mean(vdem2, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, 
                         grp = cown)

vdem3 <- sjmisc::de_mean(vdem3, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, v2xps_party,  
                         grp = cown)


#### Model Building: Complex REWB Models ####

table(vdem$gwf_type_mp)

## Relevel Factor's ##
vdem$gwf_type <- as.factor(vdem$gwf_type)
vdem<- vdem %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))
vdem$gwf_type_mp <- as.factor(vdem$gwf_type_mp)
vdem<- vdem %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))


## Relevel Factor's ##
vdem2$gwf_type <- as.factor(vdem2$gwf_type)
vdem2<- vdem2 %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))
vdem2$gwf_type_mp <- as.factor(vdem2$gwf_type_mp)
vdem2<- vdem2 %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))

## Relevel Factor's ##
vdem3$gwf_type <- as.factor(vdem3$gwf_type)
vdem3<- vdem3 %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))
vdem3$gwf_type_mp <- as.factor(vdem3$gwf_type_mp)
vdem3<- vdem3 %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))


#### Main Models, M1 - M4 ####

m1<- lmer(gender_inclusion_input_mean ~ year + gwf_type + 
            (1 + year | country_name), 
          data = vdem, 
          REML = TRUE)
sjstats::icc(m1)

m2 <- lmer(gender_inclusion_input_mean ~ year + gwf_type_mp + 
             (1 + year | country_name), 
           data = vdem, 
           REML = TRUE)

m3<- lmer(gender_inclusion_input_mean ~ year + gwf_type + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1  + year | country_name), 
          data = vdem2, 
          REML = TRUE)

m4<- lmer(gender_inclusion_input_mean ~ year + gwf_type_mp + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1 + year | country_name), 
          data = vdem2, 
          REML = TRUE)
isSingular(m4, tol = 1e-05)

#### Extracting findings and plotting effects ####

tab_model(m1, m2, m3, m4,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3", "Model 4"))

library(texreg)
texreg(list(m1, m2, m3, m4), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Communist Party Regime",
                             "Monarchy",
                             "Party Regime",
                             "Personal",
                             "Communist Regime",
                             "Communist Regime with MC", 
                             "Military Regime with MC", 
                             "Monarchy", 
                             "Monarchy with MC", 
                             "Party Regime", 
                             "Party Regime with MC", 
                             "Personal Regime", 
                             "Personal Regime with MC", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:Input1", 
       longtable = TRUE)

#### Dot-Whisker Plot ####

library(dotwhisker)
library(broom.mixed)


m1_df <- tidy(m1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m2_df <- tidy(m2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m3_df <- tidy(m3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m4_df <- tidy(m4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

m1_df <- m1_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal")) %>%
  filter(term!= "Intercept" & term!="Year")

m2_df <- m2_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC")) %>%
  filter(term!= "Intercept" & term!="Year")

m3_df <- m3_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal", 
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Monarchy" | term == "Party" | term =="Personal")

m4_df <- m4_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC", 
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Communist Party with MC" | term == "Military Regime with MC" | 
           term =="Monarchy" | term == "Monarchy with MC" | term=="Party Regime" | term == "Party Regime with MC" | 
           term =="Personal Regime" | term == "Personal Regime with MC")

m_models <- rbind(m1_df, m3_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.25),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 1", "Model 3"),
                      labels = c("Model 1", "Model 3")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 1", "Model 3"),
                         labels = c("Model 1", "Model 3"))
} 
ggsave("calculations/Output/Figure3_Input.pdf", height = 14, width = 24, units= c("cm"))

m_models <- rbind(m2_df, m4_df)
plot_dot_whisker <- {dwplot(m_models,
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr()  + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.85, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 2", "Model 4"),
                      labels = c("Model 2", "Model 4")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 2", "Model 4"),
                         labels = c("Model 2", "Model 4"))
} 
ggsave("calculations/Output/Figure4_Input.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m1)
car::vif(m2)
car::vif(m3)
car::vif(m4)



####  Models with Party Institutionalization as Indepedent Variable ####

## Relevel Factors ##
vdem3$gwf_type <- as.factor(vdem3$gwf_type)
vdem3 <- vdem3 %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))

vdem3$gwf_type_mp <- as.factor(vdem3$gwf_type_mp)
vdem3 <- vdem3 %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))


m5<- lmer(gender_inclusion_input_mean ~ year + gwf_type + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)
isSingular(m5, tol = 1e-05)

m6<- lmer(gender_inclusion_input_mean ~ year + gwf_type_mp + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)

m7 <- lmer(gender_inclusion_input_mean ~ year + gwf_type + v2xps_party_dm + v2xps_party_gm + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
             e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
             e_total_resources_income_pc_gm + 
             (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
           data = vdem3, 
           REML = TRUE)

m8<- lmer(gender_inclusion_input_mean ~ year + gwf_type_mp + v2xps_party_dm + v2xps_party_gm + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)

#### Texreg Findings ####

tab_model(m5, m6, m7, m8, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 5", "Model 6", "Model 7", "Model 8"))

texreg(list(m5, m6, m7, m8), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Communist Party Regime",
                             "Monarchy",
                             "Party Regime",
                             "Personal",
                             "Party Institutionalization", 
                             "Party Institutionalization (b)", 
                             "Communist Regime",
                             "Communist Regime with MC", 
                             "Military Regime with MC", 
                             "Monarchy", 
                             "Monarchy with MC", 
                             "Party Regime", 
                             "Party Regime with MC", 
                             "Personal Regime", 
                             "Personal Regime with MC", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:Input2", 
       longtable = TRUE)


#### Dot-Whisker Plots ####

m5_df <- tidy(m5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m6_df <- tidy(m6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m7_df <- tidy(m7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m8_df <- tidy(m8, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

m5_df <- m5_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)")) %>%
  filter(term!= "Intercept" & term!="Year")

m6_df <- m6_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)")) %>%
  filter(term!= "Intercept" & term!="Year")

m7_df <- m7_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Monarchy" | term == "Party" | term =="Personal" |
           term == "Party Institutionalization" | term == "Party Institutionalization (b)")

m8_df <- m8_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 8") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC",
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Communist Party with MC" | term == "Military Regime with MC" | 
           term =="Monarchy" | term == "Monarchy with MC" | term=="Party Regime" | term == "Party Regime with MC" | 
           term =="Personal Regime" | term == "Personal Regime with MC" | term == "Party Institutionalization" | 
           term == "Party Institutionalization (b)")

m_models <- rbind(m5_df, m7_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 5", "Model 7"),
                      labels = c("Model 5", "Model 7")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 5", "Model 7"),
                         labels = c("Model 5", "Model 7"))
} 
ggsave("calculations/Output/Figure5_Input.pdf", height = 14, width = 24, units= c("cm"))

m_models <- rbind(m6_df, m8_df)
plot_dot_whisker <- {dwplot(m_models,
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr()  + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 6", "Model 8"),
                      labels = c("Model 6", "Model 8")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 6", "Model 8"),
                         labels = c("Model 6", "Model 8"))
} 
ggsave("calculations/Output/Figure6_Input.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m5)
car::vif(m6)
car::vif(m7)
car::vif(m8)

##############################################################################################
################################################################################################
################################################################################################

#### Output Side ####


vdem <- readRDS("calculations/data/vdem_merged.rds")

#### Prepare Data for analysis ####
library(scales)

summary(vdem$gender_inclusion_input_mean)
vdem$gender_inclusion_input_mean <- rescale(vdem$gender_inclusion_input_mean, to = c(0,1))
summary(vdem$gender_inclusion_input_mean)

summary(vdem$gender_inclusion_output_mean)
vdem$gender_inclusion_output_mean <- rescale(vdem$gender_inclusion_output_mean, to = c(0,1))
summary(vdem$gender_inclusion_output_mean)

# delete those countries with less than 3 country-year observations to dela with singular fit 

vdem <- vdem %>%
  drop_na(cown, year, gwf_type, gender_inclusion_output_mean)

vdem <- vdem %>% # delete those countries with less than 3 country-year observations to dela with singular fit 
  dplyr::group_by(cown) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(cown)

vdem2 <- vdem %>%
  drop_na(cown, year, gwf_type, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio)

vdem3 <- vdem %>%
  drop_na(cown, year, gwf_type, e_wb_pop_ln, e_migdppcln, e_total_resources_income_pc, eprratio, v2xps_party)

#### Computing group_mean and de-meaned variables ####

vdem2 <- sjmisc::de_mean(vdem2, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, 
                         grp = cown)

vdem3 <- sjmisc::de_mean(vdem3, e_wb_pop_ln, e_migdppcln, eprratio, e_total_resources_income_pc, v2xps_party,  
                         grp = cown)


#### Model Building: Complex REWB Models ####

table(vdem$gwf_type_mp)

## Relevel Factor's ##
vdem$gwf_type <- as.factor(vdem$gwf_type)
vdem<- vdem %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))
vdem$gwf_type_mp <- as.factor(vdem$gwf_type_mp)
vdem<- vdem %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))


## Relevel Factor's ##
vdem2$gwf_type <- as.factor(vdem2$gwf_type)
vdem2<- vdem2 %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))
vdem2$gwf_type_mp <- as.factor(vdem2$gwf_type_mp)
vdem2<- vdem2 %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))

## Relevel Factor's ##
vdem3$gwf_type <- as.factor(vdem3$gwf_type)
vdem3<- vdem3 %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))
vdem3$gwf_type_mp <- as.factor(vdem3$gwf_type_mp)
vdem3<- vdem3 %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))

#### Main Models, M1 - M4 ####

m1<- lmer(gender_inclusion_output_mean ~ year + gwf_type + 
            (1 + year | country_name), 
          data = vdem, 
          REML = TRUE)
sjstats::icc(m1)

m2 <- lmer(gender_inclusion_output_mean ~ year + gwf_type_mp + 
             (1 + year | country_name), 
           data = vdem, 
           REML = TRUE)

m3<- lmer(gender_inclusion_output_mean ~ year + gwf_type + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1  + year | country_name), 
          data = vdem2, 
          REML = TRUE)

m4<- lmer(gender_inclusion_output_mean ~ year + gwf_type_mp + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1 + year | country_name), 
          data = vdem2, 
          REML = TRUE)

#### Extracting findings and plotting effects ####

tab_model(m1, m2, m3, m4,
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 1", "Model 2", "Model 3", "Model 4"))

library(texreg)
texreg(list(m1, m2, m3, m4), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Communist Party Regime",
                             "Monarchy",
                             "Party Regime",
                             "Personal",
                             "Communist Regime",
                             "Communist Regime with MC", 
                             "Military Regime with MC", 
                             "Monarchy", 
                             "Monarchy with MC", 
                             "Party Regime", 
                             "Party Regime with MC", 
                             "Personal Regime", 
                             "Personal Regime with MC", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:Output1", 
       longtable = TRUE)

#### Dot-Whisker Plot ####

library(dotwhisker)
library(broom.mixed)


m1_df <- tidy(m1, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m2_df <- tidy(m2, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m3_df <- tidy(m3, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m4_df <- tidy(m4, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

m1_df <- m1_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 1") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal")) %>%
  filter(term!= "Intercept" & term!="Year")

m2_df <- m2_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 2") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC")) %>%
  filter(term!= "Intercept" & term!="Year")

m3_df <- m3_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 3") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal", 
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Monarchy" | term == "Party" | term =="Personal")

m4_df <- m4_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 4") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC", 
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Communist Party with MC" | term == "Military Regime with MC" | 
           term =="Monarchy" | term == "Monarchy with MC" | term=="Party Regime" | term == "Party Regime with MC" | 
           term =="Personal Regime" | term == "Personal Regime with MC")

m_models <- rbind(m1_df, m3_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.25),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 1", "Model 3"),
                      labels = c("Model 1", "Model 3")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 1", "Model 3"),
                         labels = c("Model 1", "Model 3"))
} 
ggsave("calculations/Output/Figure3_Output.pdf", height = 14, width = 24, units= c("cm"))

m_models <- rbind(m2_df, m4_df)
plot_dot_whisker <- {dwplot(m_models,
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr()  + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.85, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 2", "Model 4"),
                      labels = c("Model 2", "Model 4")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 2", "Model 4"),
                         labels = c("Model 2", "Model 4"))
} 
ggsave("calculations/Output/Figure4_Output.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m1)
car::vif(m2)
car::vif(m3)
car::vif(m4)



####  Models with Party Institutionalization as Indepedent Variable ####

## Relevel Factors ##
vdem3$gwf_type <- as.factor(vdem3$gwf_type)
vdem3 <- vdem3 %>% 
  mutate(gwf_type = relevel(gwf_type, "military"))

vdem3$gwf_type_mp <- as.factor(vdem3$gwf_type_mp)
vdem3 <- vdem3 %>% 
  mutate(gwf_type_mp = relevel(gwf_type_mp, "military regime"))


m5<- lmer(gender_inclusion_output_mean ~ year + gwf_type + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)
isSingular(m5, tol = 1e-05)

m6<- lmer(gender_inclusion_output_mean ~ year + gwf_type_mp + v2xps_party_dm + v2xps_party_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)

m7 <- lmer(gender_inclusion_output_mean ~ year + gwf_type + v2xps_party_dm + v2xps_party_gm + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
             e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
             e_total_resources_income_pc_gm + 
             (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
           data = vdem3, 
           REML = TRUE)

m8<- lmer(gender_inclusion_output_mean ~ year + gwf_type_mp + v2xps_party_dm + v2xps_party_gm + e_migdppcln_dm + e_migdppcln_gm + e_wb_pop_ln_dm +
            e_wb_pop_ln_gm +  eprratio_dm + eprratio_gm + e_total_resources_income_pc_dm + 
            e_total_resources_income_pc_gm + 
            (1 + year | country_name) + (1 + v2xps_party_dm | country_name), 
          data = vdem3, 
          REML = TRUE)

#### Texreg Findings ####

tab_model(m5, m6, m7, m8, 
          show.ci = FALSE, 
          show.se = TRUE, 
          auto.label = FALSE, 
          string.se = "SE",
          show.icc = FALSE, 
          dv.labels = c("Model 5", "Model 6", "Model 7", "Model 8"))

texreg(list(m5, m6, m7, m8), 
       head.tag = TRUE, body.tag = TRUE,
       digits = 2,
       custom.coef.names = c("(Intercept)",
                             "Year",
                             "Communist Party Regime",
                             "Monarchy",
                             "Party Regime",
                             "Personal",
                             "Party Institutionalization", 
                             "Party Institutionalization (b)", 
                             "Communist Regime",
                             "Communist Regime with MC", 
                             "Military Regime with MC", 
                             "Monarchy", 
                             "Monarchy with MC", 
                             "Party Regime", 
                             "Party Regime with MC", 
                             "Personal Regime", 
                             "Personal Regime with MC", 
                             "GDP pc log",
                             "GDP pc log (b)",
                             "Population log",
                             "Population log (b)",
                             "Ethnic Excluded Population", 
                             "Ethnic Excluded Population (b)", 
                             "Resource Income",
                             "Resource Income (b)"),
       booktabs = TRUE,
       use.packages = FALSE,
       caption = "Linear Within-Between Model Predicting Women's Political Exclusion",
       fontsize = "scriptsize",
       stars = c(0.001, 0.01, 0.05, 0.1), 
       symbol = "\\dagger",
       label = "Table:Output2", 
       longtable = TRUE)


#### Dot-Whisker Plots ####

m5_df <- tidy(m5, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m6_df <- tidy(m6, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m7_df <- tidy(m7, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")
m8_df <- tidy(m8, conf.int=TRUE, exponentiate=FALSE, effects="fixed",
              conf.method = "Wald")

m5_df <- m5_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 5") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)")) %>%
  filter(term!= "Intercept" & term!="Year")

m6_df <- m6_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 6") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)")) %>%
  filter(term!= "Intercept" & term!="Year")

m7_df <- m7_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 7") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_typecommunist party" = "Communist Party", 
                       gwf_typemonarchy = "Monarchy", 
                       gwf_typeparty = "Party", 
                       gwf_typepersonal = "Personal", 
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Monarchy" | term == "Party" | term =="Personal" |
           term == "Party Institutionalization" | term == "Party Institutionalization (b)")

m8_df <- m8_df %>%
  by_2sd(vdem2) %>%
  arrange(term) %>%
  mutate(model="Model 8") %>%
  relabel_predictors(c("(Intercept)" = "Intercept",                       
                       year = "Year", 
                       "gwf_type_mpcommunist party" = "Communist Party", 
                       "gwf_type_mpcommunist party with ME" = "Communist Party with MC", 
                       "gwf_type_mpmilitary regime with ME" = "Military Regime with MC", 
                       "gwf_type_mpmonarchy" = "Monarchy", 
                       "gwf_type_mpmonarchy with ME" = "Monarchy with MC",
                       "gwf_type_mpparty regime" = "Party Regime", 
                       "gwf_type_mpparty regime with ME" = "Party Regime with MC",
                       "gwf_type_mppersonal regime" = "Personal Regime", 
                       "gwf_type_mppersonal regime with ME"= "Personal Regime with MC",
                       v2xps_party_dm = "Party Institutionalization", 
                       v2xps_party_gm = "Party Institutionalization (b)",
                       e_migdppcln_dm =  "GDP pc log (within)", 
                       e_migdppcln_gm = "GDP pc log (between)", 
                       e_wb_pop_ln_dm = "Population (within)", 
                       e_wb_pop_ln_gm = "Population (between)", 
                       eprratio_dm = "Excluded Population (within)", 
                       eprratio_gm = "Excluded Population (between)", 
                       e_total_resources_income_pc_dm = "Resource income (within)", 
                       e_total_resources_income_pc_gm = "Resource income (between)")) %>%
  filter(term == "Communist Party" | term=="Communist Party with MC" | term == "Military Regime with MC" | 
           term =="Monarchy" | term == "Monarchy with MC" | term=="Party Regime" | term == "Party Regime with MC" | 
           term =="Personal Regime" | term == "Personal Regime with MC" | term == "Party Institutionalization" | 
           term == "Party Institutionalization (b)")

m_models <- rbind(m5_df, m7_df)
plot_dot_whisker <- {dwplot(m_models, 
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr() + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 5", "Model 7"),
                      labels = c("Model 5", "Model 7")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 5", "Model 7"),
                         labels = c("Model 5", "Model 7"))
} 
ggsave("calculations/Output/Figure5_Output.pdf", height = 14, width = 24, units= c("cm"))

m_models <- rbind(m6_df, m8_df)
plot_dot_whisker <- {dwplot(m_models,
                            dot_args = list(size = 3, aes(shape = model)),
                            whisker_args = list(size = 1, aes(linetype = model)),
                            dodge_size = 0.5) +
    theme_pubr()  + xlab("Coefficient Estimate") + ylab("") +
    geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
    ggtitle("") +
    theme(plot.title = element_text(face="bold"),
          legend.position = c(0.05, 0.8),
          legend.justification = c(0, 0), 
          legend.background = element_rect(colour="grey80"),
          legend.title = element_blank()) +
    scale_colour_grey(start = .4, end = .1, # if start and end same value, use same colour for all models 
                      name = "Model", 
                      breaks = c("Model 6", "Model 8"),
                      labels = c("Model 6", "Model 8")) +
    scale_shape_discrete(name = "Model",
                         breaks = c("Model 6", "Model 8"),
                         labels = c("Model 6", "Model 8"))
} 
ggsave("calculations/Output/Figure6_Output.pdf", height = 14, width = 24, units= c("cm"))

## Testing Variance Inflation factor for multicollinearity ##

car::vif(m5)
car::vif(m6)
car::vif(m7)
car::vif(m8)

